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SUMMARY 

We conducted numerical simulations of the wind-forcing of the sea level variations in the North Sea 
using a barotropic ocean model with realistic geography, bathymetry, and boundary conditions, to 
examine the forcing of the 14-month “pole tide” which is known to be strong along the Denmark- 
Netherlands coast. The simulation input is the monthly-mean surface wind stress field from the 
National Centers for Environmental Prediction (NCEP) reanalysis for the 40-year period 1958- 1997 . 
The output sea level response was then compared with 10 coastal tide gauge records from the 
Permanent Service for Mean Sea Level (PSMSL). Besides the strong seasonal variations, several 
prominent quasi-periodicities exist at around 7 years, 3 years, 14 months, 9 months, and 6.5 months. 
Correlation and spectral analyses show remarkable agreement between the model output and the 
observations, particularly in the 14-month, or Chandler period band. The latter indicates that the 
enhanced pole tide found in the North Sea along the Denmark-Netherlands coast is actually the 
coastal setup response to wind stress forcing with a periodicity of 14 months. We find no need to 
invoke a geophysical explanation involving resonance-enhancement of pole tide in the North Sea to 
explain the observations. 

1. INTRODUCTION 

In this paper we reconsider the following question: Is the anomalously large 14-month pole tide 
observed at the tide gauge stations along the continental coast of the North Sea really the pole tide 
that is somehow enhanced, or is it primarily of direct meteorological origin (and hence not really the 
“pole tide”) forced by the regional surface wind stress with a 14-month periodicity? Our results 
strongly support the latter. 

The polar motion of the earth’s rotational axis sets up the pole tide in the ocean via the centrifugal 
potential induced by the polar motion. In particular, the 14-month Chandler wobble in the polar 
motion generates a 14-month pole tide. At this long period, there is little horizontal motion associated 
with the pole tide and the height is close to the equilibrium value, or height the tide would have if all 
the wobble poten tial energy were used to raise the sea level. Theoretical studies by Dickman (1985, 
1988) and others have shown that the pole tide in the deep oceans should be near its equilibrium 
value. The maximum amplitude of an equilibrium pole tide is about 0.5 cm at 45° latitude. 
Independent data studies over the past several decades (for a review, see, e.g., Lambeck 1980, p. 
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212-214) have shown that while the pole tide amplitude is generally near the equilibrium tide value, 
it is anomalously large, up to 3 cm, in parts of the North and Baltic Seas. A salient feature in the 
North Sea is the enhanced 14-month tide along the continental coast known as the eastward 
intensification (e.g., Xie & Dickman 1995), where the tidal amplitude increases sharply eastward 
along the Denmark-Netherlands coast to five times its equilibrium value (Miller & Wunsch 1973, Fig. 
2, reproduced in Lambeck, 1980, Fig. 8.1). A number of studies, most recently by Ekman & 
Stigebrandt (1990), have concluded that the secular change in the pole tide, especially in the North 
and Baltic Seas, is not correlated with the pole tide forcing, and have suggested a meteorological 
component to the forcing. 

Analytical and numerical studies of the dynamics of the pole tide in the North Sea have been made 
by Wunsch (1974, 1986), Dickman & Preisig (1986), Carton & Wahr (1986), and Xie & Dickman 
( 1995). These are mentioned by Wunsch et al. (1997) who discuss the possibility of a basin resonance 
in the North Sea at the pole tide period. No satisfactory explanation has yet been given for the 
enhancement of the North Sea pole tide along the continental coast as a result of the direct pole tide 
forcing, the forcing by the global tide as an open boundary condition, or any possible resonance 
phenomenon with a 14-month period. These possibilities were considered and rejected years ago by 
Miller (1973), who then suggested that the wind stress was the most probable cause of the observed 
14-month tide. Indeed, these studies collectively show that a pole tide enhancement cannot be derived 
from any homogeneous system of equations describing the flow regime in the North Sea. An analytic 
study by O'Connor (1986) showed that the eastward intensification can be explained in terms of an 
inhomogeneous (forced) system of equations with a linear dynamical model for an idealized 
rectangular basin, assuming that the forcing is a small change in the east-west wind stress component 
near the 14-month Chandler period. 

The North Sea climate is strongly influenced by the changes in the North Atlantic wind field 
(WASA Group 1998). There is evidence for a periodic 14-month oscillation in sea surface pressure 
over the North Atlantic, North Sea, and Baltic Sea regions, which is directly related to an oscillation 
in the position of the Icelandic low as part of the North Atlantic Oscillation system (e.g., Wallace & 
Gutzler 1981). This was demonstrated by Bryson & Starr (1977) who analyzed a National Center for 
Atmospheric Research (NCAR) global monthly mean sea level pressure data set for the period 
1900-1970. Other previous investigations into this phenomenon are listed by O'Connor (1986) and 
Tsimplis, Flather & Vassie (1994). Two observational studies investigated periodicities in the 
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14-month meteorological forcing over the North Sea region and their correlation with the pole tide 
observations there. The first study was by Trupin & Wahr (1990) who used the NCAR monthly mean 
sea surface pressure data. They concluded that the pole tide in the North Sea was correlated with 
the atmospheric pressure at that 14-month period. It might have been more appropriate to make a 
correlation with the north-south pressure gradient in that region, since this determines the east-west 
component of the wind through the geostrophic relation, modified by friction in the boundary layer. 
In the second study, Tsimplis et al. (1994) used a 30-year time series of wind and pressure data from 
the Norwegian Meteorological Institute and tidal specifications on the northern open boundary to run 
a numerical tide and surge model of the North Sea. They concluded that the pole tide along the 
continental coast was the result of a variation in the east- west wind stress with a 14-month 
periodicity. However, there remains some controversy on whether the amplitude and phase of the 
wind forcing spectral peak over the North Sea is sufficient to account for the sea level response near 
the Chandler period. This is reflected in the comments of Wunsch et al. (1997) on the work of 
Tsimplis et al. (1994). 

We shall re-investigate this problem with a new, and hopefully more complete meteorological data 
set and a North Sea ocean model. Historically, meteorological data sets have been compiled only 
from observations. More recently, output of global numerical weather prediction models assimilating 
all available observational data are being used to make climatologies. The data for the surface winds 
are accepted as generally accurate because they are representative of the large scale flow. We shall 
adopt the 40-year reanalysis product of the National Centers for Environmental Prediction (NCEP) 
(Kalnay et al. 1996). Figure 1 depicts the geography of the region studied, as well as its average 
NCEP surface wind stress field (see below), which we use to force a barotropic North Sea ocean 
model generally following the approach of Tsimplis et al. (1994). We compare our model output, or 
“computed” sea level with actual observations for 10 continental coastal stations from the Permanent 
Service for Mean Sea Level (PSMSL). Good agreement is found, including in the Chandler frequency 
band, indicating that the enhanced “pole tide” in the North Sea is specifically the coastal setup 
response to wind stress forcing with a periodicity of 14 months. 

2. OBSERVED SEA LEVELS 

The PSMSL (URL http://www.nbi.ac.uk/psmsl/psmsl.info.html) archives monthly sea level records 
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from tide gauge stations worldwide. For the North Sea, it provides 17 station records that started 
prior to 1958 (and many more records that are shorter). Thirteen of them are on the continental coast, 
and we examine 10 of these stations (Fig. 1) from the English Channel to Denmark, where the pole 
tide amplitude becomes pronounced. They are, from north to south: Esbjerg, West-Terschelling, 
Delfzijl, Harlingen, Den Helder, Ijmuiden, Hoek Van Holland, Maassluis, Vlissingen, Oostende (see 
Table 1). The northern stations are known to have the strongest 14-month tide showing the eastward 
intensification; the first five stations are thus combined later for the purpose of studying this signal. 

3. COMPUTED SEA LEVELS 

3. 1 Atmospheric wind stress forcing 

A description of the monthly mean wind directions and persistency at various locations over the 
North Sea is given by Korevaar (1990). During the period 1961-1980 the prevailing annual wind 
direction was from the SW to WNW. The seasonal variation was such that in summer the wind was 
from WNW to NW, while in winter it was from SW to W. In the spring and fall, the winds can be 
from several directions, and the monthly wind direction can have an easterly component. Also, the 
strongest winds usually have a westerly component. This information is well known, but here we 
want to make use of the latest meteorological data to search for interested periodicities in the wind 
forcing of the North Sea. 

The global data set of the NCEP reanalysis (Kalnay et al. 1996) for the 40-year period Jan. 1958 - 
Dec. 1997 contains monthly mean data for a large number of variables. We first examine directly the 
low level 1000 millibar (mb) winds in the air-sea boundary layer. The NCEP data set has the monthly 
mean values for the eastward (u) and northward (v) horizontal components, given at 2.5° x 2.5° 
intervals. Those points over the North Sea and immediately surrounding region that are used in this 
study are at the intersection of longitudes 5.0° W, 2.5° W, 0°, 2.5° E, 5.0° E, 7.5° E, 10.0° E, and 
latitudes 50.0° N, 52.5° N, 55.0° N, 57.5° N, 60.0° N. 

Figure 2 summarizes the behavior of the NCEP u and v wind fields of the North Sea region, 
averaged over the above grid points. Fig. 2(b) gives their maximum-entropy power spectra, where 
the order of the corresponding autoregressive model is selected according to criteria in Marple 
(1990), about 20% of the record length in our case. While not surprisingly both show strong annual 
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power, the u component has a rich ensemble of periodicities, including a prominent signal at the 
Chandler period near 14 months. The latter, as we show later, is the major driving source of the 
observed 14-month signal in the sea leveL The v component did not exhibit strong periodicities other 

than the annual. 

The parameter we use to directly force the ocean model is the horizontal wind stress vector on 
the sea surface, t = (t x , x y ), which is a quadratic function of the surface wind velocity V=(u,v). 

T= p C d V IV! (1) 


where p is the surface layer air density and C d is the drag coefficient. The air density varies only 
sli ghtly with air temperature and humidity. The drag coefficient is a function of both the wind speed 
and the atmospheric surface layer stability. The air-sea temperature difference and heat flux increase 
momentum transfer due to turbulence at the air-sea interface. These stability effects can cause the 
drag coefficient to significantly vary seasonally and spatially. The NCEP global model uses a surface 
layer parameterization to calculate the wind stress (e.g. Haltiner & Williams 1980, section 8-6). The 
drag coefficient is calculated by the formula C d = k 2 / [ln(z/z c ) - t|r(z/L)] 2 , where k=0.4 is the von 
Karman constant, z is the first model level (usually 30 - 40 m), z c is the roughness length, L is the 
Monin-Obukhov length scale which depends on turbulent fluxes, and t|r is an empirical stability 
function. This latter depends on the Richardson number, which relates the turbulence production due 
to wind shear to the turbulence damping due to the stratification, resulting from air-sea temperature 
differences. In fact, the value of the 10 m wind velocity can be found from V = (V*/k) [ln(z/z 0 ) - 
t|i(z/L)], where V* is the friction velocity, which can be related to the roughness length over water 
by the Chamock relation z Q = 0.014 V* 7g, where g is the acceleration of gravity. The equations are 


solved numerically by iteration. 

The NCEP wind stress data are given at intervals of 1.875° in longitude and 1.904° in latitude. 


The points over the North Sea and surrounding region that are used in this study are at the 
intersections of longitudes 3.75° W, 1.875° W, 0°, 1.875° E, 3.75° E, 5.625 E, 7.5 E, 9.375 E, 
and latitudes 48.571° N, 50.475° N, 52.38° N, 54.285° N, 56.189° N, 58.094“ N, 59.99° N. Fig. 1 


shows the “prevailing”, 40-year average of these monthly NCEP wind stress fields at these grid 
points. This compares favorably with the annual mean wind stress over the North Sea given by 
Hellerman (1967), who found some spatial variability in the annual wind stress pattern with the 
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representative values t x = 0.062 nt m' 2 , t y = 0.021 nt m' 2 . 

The NCEP wind stresses were calculated at 40-minute time step intervals, from which a daily 
average was made. The daily wind stress output of the NCEP model was then averaged to form the 
monthly mean values. Note that the monthly averaging must be invoked on the wind stress values 
after they are computed from the wind velocity, rather than on the wind velocity itself prior to 
calculation of wind stress. This is because the averaging of wind velocity would tend to reduce the 
resultant wind stress due to the quadratic nature of Eq. 1 (see below for more discussion). 

We made an alternate wind stress calculation as a check on the NCEP surface wind stress output. 
This was for neutral stability conditions, where the air-sea temperature difference was neglected. This 
has been done in previous studies when the air-sea temperatures and stability conditions were 
unobtainable. The wind stress was calculated directly from equation (1) with the NCEP 1000 mb 
winds used in place of 10 m surface winds, a neutral drag coefficient, and the constant air density 1.2 
kg m 3 . A climatology of the global ocean wind stress based on ECMWF 1000 mb winds was 
produced by Trenberth, Large & Olson (1989) for the years 1980-1986, during which only the 1000 
mb winds were available. As explained by Trenberth et al. (1989), the ship wind observations were 
assigned to the 1000 mb level and effectively represent 10 m winds. However, an evaluation of 
ECMWF-based climatological wind stress fields by Mestas-Nunez et al. (1994) concluded that the 
1000 mb winds were more representative of winds at a height of 30 m, and the resulting wind stresses 
of Trenberth et al. (1989) were too large by 50%. 

The drag coefficient depends on the sea state (waves, swell) and is written as a linear function 
of wind speed. A representative value for the North Sea determined from, e.g.. Smith et al. (1992) 
is C N = (0.50 + 0.091 1 VI ) x 10 \ where the wind speed is in ms' 1 . The wind stresses calculated under 
neutral conditions with the 1000 mb mean monthly winds were smaller than the wind stresses 
produced directly by the NCEP model by about a factor of five. This is explained by the fact that 
the wind stress is proportional to the square (or even cube) of the wind speed, and so the momentum 
input is strongly influenced by the high wind events. The monthly averaging done on the 1000 mb 
wind velocity itself smooths out these high wind events (see above), and this is at least partly why our 
neutral stability calculation significantly underpredicts the wind stress. Other researchers (e.g.. 
Stammer et al. 1996) have noted that if an ocean model does not resolve eddies and have the correct 
heat flux surface conditions, it can underpredict the sea level and its variability by a factor of over 2. 

In our present application, we conclude that the NCEP wind stress calculation should be more 
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accurate. It will be adopted for our ocean model forcing below. 


3.2 North Sea dynamics and model 


The flow in the North Sea has been studied using numerical models (see Sundermann & Lenz 1983 
for a review). These models have had great success in modeling the sea level and currents in response 
to tidal and wind forcing. They include both realistic coastlines and bottom topography. Tides and 
wind forced surges in the North Sea can be modeled adequately with barotropic models, which solve 
for the sea level and the vertically averaged currents (Flather 1984; Verboom, de Ronde & van Dikj 
1992). The numerical model used by Tsimplis et al. (1994) is described by Flather, Proctor & Wolf 
(1991). We will use the numerical ocean model developed at Princeton University and described by 
Blumberg & Mellor (1987), written in generalized coordinates. Since we desire to investigate the 
wind forced sea level set up along the continental coast, we can use the barotropic subset of these 
model equations. These are the flux form of the vertically integrated nonlinear continuity and 
horizontal momentum equations, which are solved for the sea level and the horizontal components 
of velocity. In spherical coordinates, the horizontal momentum equations consist of terms for the 
local time derivative and horizontal advective terms, the Coriolis acceleration, sea level pressure 
gradient, tangential wind stress on the sea surface, and quadratic bottom friction with a constant drag 
coefficient of 0.0025. The system of equations is solved using a finite differencing scheme that is 

centered in time and space on the Arakawa-C grid. 

The model North Sea domain (Fig. 1) covers the region from 51° N to 59° N latitude, and from 
3° W to 9° E longitude. There are 97 grid points in the north-south with a spacing of 5' in latitude 
(9.27 km), and 73 grid points in the east- west with a spacing of 10' in longitude (9.55 - 1 1.66 km). 
The bathymetry was digitized from the British Admiralty charts 2 182 A, 2182B, and 2182C. Over 
most of the North Sea the depths were taken at intervals of 15' latitude by 30' longitude and then 
interpolated to the finer grid with a bilinear interpolation scheme. However, in order to accurately 
model the wind forced setup along the European coast, the depths either south of 54° N latitude or 
east of 7° E longitude were digitized at intervals of 5' latitude and 10 longitude. The grid size did 
not permit resolution of the barrier islands and channels to ports along the continental coast. Thus 
comparisons with tide gage data at some ports were made with the water level data from the closest 
model grid points. The shallowest depths recorded were 3 m along the continental coast, while the 
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deepest was 570 m in the Norwegian Trench. Any grid point depths shallower than 100 m along the 
Norwegian coast were made land points. The domain was closed along the southern, western and 
eastern boundaries to avoid problems with unnecessary open boundaries. A time step of 60 s is short 
enough to ensure computational stability. 

We now describe the dynamics of a residual circulation, which is the mean velocity field averaged 
over a period of at least several days, to cancel transitory wind currents and tidal oscillations. The 
flow in the North Sea can be influenced by the open boundary conditions, the atmospheric pressure 
gradient (inverted barometer effect), and the wind stress. The work of Timmerman (1977, p.47-49), 
Prandle (1975), and Flather (1984) indicates that the sea level and circulation in the southern North 
Sea are predominantly the result of the wind stress, rather than the atmospheric pressure gradient. 
Tsimplis et al. (1994) included the semi-diurnal and diurnal tidal forcing as open boundary conditions 
and showed that this had no effect on the 14-month tide. Carton & Wahr (1986) showed that using 
the deep ocean pole tide as an open boundary condition resulted in raising the elevation of the basin 
without a sea level gradient along the continental coast. Accordingly, only the wind stress forcing 
will be considered in this study. 

When the winds are steady or slowly varying in time, the model is not sensitive to the open 
northern boundary conditions which is removed from the region of our main interest, the continental 
Denmark-Netherlands coast. The equations are linearized at the northern boundary, and a radiation 
condition is used there: v = (g/H) ,/2 r) - t x / ( p 0 f H ), where H is the depth, r| is the sea level, p 0 
is the seawater density, and f is the value of the Coriolis parameter at 59° N latitude. We have seen 
that the seasonal and 14- month periodic wind stresses over the North Sea have a strong eastward 
component. This boundary condition is particularly well suited for modeling the flow due to winds 
with a strong eastward component, since it allows for southward Ekman drift into the North Sea 
across most of the northern boundary, and northward geostrophic outflow near the eastern boundary, 
as discussed by O'Connor (1986). 

3.3 Model tests 

The North Sea model used is a one layer barotropic hydrodynamic model forced only by the 
wind stress and does not consider density changes. The amplitude of the ocean model sea level 
depends on the model parameters (the bathymetry and coastline, bottom friction coefficient, and open 
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boundary conditions at the northern boundary), in addition to the wind stress forcing. The model grid 
does not resolve the barrier islands and channels to ports along the continental coast. In some cases 
the tide gauges are at ports inside of barrier islands. Presumably, the amplitudes could be improved 
by engineering the drag coefficient, friction, and station depths until better agreement is reached, as 
is typically done for storm surge models. In the ocean the annual cycle in sea level is the result of 
both the wind stress forcing and the steric effects due to density changes from temperature and 
salinity variations. 

Despite this limitation of the barotropic model, one way to tell how "good" it is can be seen by 
running the model with a constant annual average wind stress, and looking at the setup in sea level 
akmg the continental coast. As a first test, the North Sea model was run with this constant 40-year 
average NCEP wind stress forcing (Fig. 1). The average stress for NCEP grid points were 
interpolated to the ocean model grid with a bilinear interpolation scheme. The model was then run 
until a steady state resulted after three days. The resulting sea levels and currents are shown in Fig. 
3, with an evident eastward intensification of sea level along the continental coast. In particular there 
is a maximum setup of about 12 cm along the Denmark coast. These results are in very close 
agreement with the results of Prandle (1978) and Davies (1982), who forced models of the North Sea 
with an average annual wind stress. These results are also in agreement with the observational studies 
of Pattullo et al. (1955) and Rossiter (1967), which show that the mean annual sea level and 
amplitude of the seasonal oscillation exhibit a steep gradient along the Netherlands coast, which they 

attributed to eastward winds and atmospheric pressure. 

The residual circulation that results from these eastward winds is a counterclockwise gyre over 
the North Sea, as a consequence of the slope of the bottom topography. This has been shown by 
analytical models with north-south bottom slope (Fumes 1980, O’Connor 1986). However, realistic 
bathymetry is necessary to obtain the strong sea level gradients along the continental coast. In the 
deeper northern part, there is inflow due to Ekman transport. Water piles up in the shallow German 
Bight and flows northward along the continental coast in a geostrophic boundary current. This gives 
rise to a sea level gradient along the Netherlands coast. This is shown clearly in the numerical 
modeling studies of Maier-Reimer (1977), Prandle (1978; 1984), Pingree & Griffiths (1980), Davies 
(1982), and Backhaus & Maier-Reimer (1983). It should be noted that storm surge models of the 
North Sea show that when there are strong persistent westerly winds for a day or longer, there is a 
counterclockwise circulation and sea level setup similar to that described above, with strong sea level 
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gradient along the continental coast (Prandle 1975, Timmerman 1977, Mathisen & Johansen 1983, 
and Bijl 1997). 

In this paper, we only consider the model response to monthly-mean wind stress. Since the 
response of the shallow North Sea to changing winds is the period of the first external mode sea level 
oscillation, on the order of several days, we may consider a mean monthly flow to be in a steady state. 
This approach is justified because we are investigating the sea level set up by wind forcing on 
seasonal and longer time scales, which vary much more slowly than the adjustment time of the North 
Sea to changes in winds. For each wind stress forcing that is used, the model is spun up from rest 
by increasing the wind stress forcing linearly from zero (ramping) over one inertial period (0.61 days) 
and then keeping it constant. When the spatially constant annual average wind stress forcing was 
used, the resulting model flow regime was almost in a steady state balance after one day, as evidenced 
by plots of currents and sea level at one, three, and five days. For this reason, it was decided that 
three days of simulation are enough to bring the model to a steady state with any monthly wind stress 
forcing. 

Thus, a three-day simulation with steady wind stress forcing was made for each of the 480 
individual months Jan. 1958 through Dec. 1997. The values of the model water levels at the grid 
points closest to selected coastal stations (see Fig. 1) were written out at the end of each three-day 
simulation to represent the value for that month. 


4. SEA LEVEL COMPARISONS 

4. 1 Seasonal signals 

We now compare our North Sea model output, or the computed sea levels with those actual 
observed at the 10 stations mentioned above. We first examine the seasonal signals at the annual and 
semi-annual periods. To estimate the amplitude and phase of the seasonal signals, we simultaneously 
least-squares fit sine and cosine functions with the annual and semi-annual periods to the 40-year time 
series. The resultant sine and cosine functions are then converted into the form 
(amplitude) *sin(cot+phase), where t=0 refers to the nominal epoch of January 15.5. 

Fig. 4 shows a comparison between the computed and observed annual and semi-annual 
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amplitude/phase as a function of latitude for the ten coastal stations. It is seen that. 

(i) The observed annual amplitudes are between 5-15 mm, increasing almost monotonically towards 
the north and east. The corresponding computed amplitudes follow the same trend but are smaller 
by about a factor of 3. 

(ii) The observed semi-annual amplitudes are between 1-3 mm, also increasing towards the north and 
east. The corresponding computed amplitudes follow the same trend but are smaller by a factor of 
about 2. 

(iii) In general the computed seasonal phase agrees well with the observed for both annual and semi- 
annual signals, except for the annual phase of the southernmost stations. 

Some of the main features above have been discussed previously by Tsimplis et al. (1994) and 
Xie & Dickman (1996). The difference in the seasonal amplitudes is expected because our computed 
values are only those caused by surface wind stress. There are other important factors that contribute 
significantly to the seasonal signals that are observed. The most important is presumably the steric 
effect due to thermal and salinity variations particularly in the surface water layers. The strength of 
air pressure effects on the sea level, such as in the form of the inverted-barometer effect, also varies 
seasonally. Furthermore, the annual and semi-annual ocean tides due to the sun s tidal force have 
centimeter level amplitudes. Another interesting possible contributor is actually the pole tide: Since 
the Earth’s polar motion contains an annual wobble besides the 14-month Chandler wobble and nearly 
equal in amplitude, the annual wobble should presumably generate a pole tide of a half centimeter 

at maximum just like the Chandler wobble would. 

It is interesting, then, that the seasonal phase estimates agree quite well, and that the computed 
and observed seasonal amplitudes show the same geographical trend. This presumably indicates that 
the sum of other seasonal effects are highly correlated with that of the wind-forcing contribution. 

4.2 Non-seasonal, broad-band signal 

Next, we remove by subtraction the seasonal signals from all (monthly) time series, in order to 
concentrate on non-seasonal, broad-band signals. We also remove any linear trend in the time series, 
whatever the cause of the trend may be. The latter may include eustatic sea level rise, post-glacial or 
local movements of the ground, or even systematic observational errors. These are outside the interest 
of the present study. 
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Comparing such non-seasonal broad-band signals, we find the following: 

(i) The computed sea level variation agrees remarkably well with the observed at the northern 
stations, while it becomes somewhat larger than the observed in overall amplitude for the southern 
stations. This is evident in Fig. 5 (and Table 1 below). 

(ii) Table 1 lists the correlation coefficients between the computed and the observed time series. The 
geographical trend is obvious: the correlation is very high for the northern stations, while decreases 
rapidly toward the southernmost stations. The cross-correlation functions (not shown) all look 
characteristically similar to Fig. 7 below. 

(iii) Since the cross-correlation only compares the “pattern” of the time series and conveys no 
information about the overall size of the amplitude, we should calculate the ratio [standard deviation 
of computed amplitude]/[standard deviation of observed amplitude]. The latter are also listed in Table 
1. The ratio is close to one, except for the few southernmost stations. 

Findings (ii) and (iii) confirm (i) in a quantitative manner. We conclude here that, while it tends 
to over-prediets the amplitude for the southernmost stations, our North Sea model using NCEP wind 
stress input predicts well the non-seasonal sea level variations along the Denmark-Netherlands coast. 
The latter is indeed where the strongest 14-month tide is found. 

Thus, we now concentrate on the five northernmost stations. We simply average the five time 
series into a “composite” series, and make comparisons between the computed and the observed sea 
levels. Fig. 6 shows the two composite series; remarkable matching is already evident. Fig. 7 gives 
their (time-domain) cross-correlation as a function of time shift between the two series. A sharp peak 
at zero time-shift represents a correlation coefficient as high as 0.86. The correlation drops to near- 
zero with time shifts of only 1 month (the time resolution of our study) and longer. A spectral 
decomposition of this correlation is given in Fig. 8. It shows the (frequency-domain) coherence 
spectrum of the two composite series. The coherence spectrum is computed using the multi-taper 
spectral technique (Thomson 1982; also, e.g., Chao & Eanes 1995), which provides robust, 
minimum-leakage spectral estimates. Seven orthogonal tapers with time-bandwidth product of 4n 
were adopted. The coherence across the entire spectrum is much higher than, say, the 99% 
confidence level (0.54 in the present case, see Chao & Eanes 1995). The phase difference remains no 
more than ±20°. 

For the same reason as in Finding (iii) above, we conduct a comparison of the individual spectra. 
Fig. 9 shows the two Fourier power spectra, again computed by means of the multi-taper technique. 
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Consistent with Tsimplis et a/.’s (1994) model results, a very good match is evident across the entire 
spectral band to the Nyquist frequency, indicating the agreement of the power, or the size of the 
amplitude of the two (composite) time series. 

Alternatively, one can show the time-varying spectral behavior of a time series by way of its time- 
frequency wavelet spectrum. The wavelet analysis was introduced to represent functions that are 
localized in time and frequency (Morlet et al. 1982). Using the Morlet wavelet, the result is shown 
in Fig. 10 (cf. Chao & Naito 1995). Here the contours indicate the strength of undulations at any 
particular frequency (or period) as a function of time going horizontally on the plot. High contrasts 
in shade indicate strong undulations; the brighter contour closures indicate positive amplitude and the 
darker ones the negative amplitude. We have to an extent saturated the grey scale in order to bring 
out the finer features. 

Let’s now examine Figs. 8-10 in more detail. In Fig. 8, a frequency band where the coherence 
dips below 99% (and the associated phase shows some “anomalous” behavior) is around the biennial 
period. This appears to be insignificant, however, as the signal power is rather low in this frequency 
band (in fact, the lowest among the entire frequency band, see Fig. 9a). It simply stems from low 
signal-to-noise ratio in this band and the resultant numerical instability. The same applies to the zero 
and near- zero frequencies, because the mean and a linear trend have been removed beforehand. 

Figs. 9 and 10 reveal some interesting periodicities in the sea level records. There are prominent 
spectral peaks and strong undulations common to both the computed and observed sea levels at 
periods of around 7 years, 3 years (especially during the first half of the records), 14 months, 9 
months, and 6.5 months. They all correspond to the peaks detected in the original eastward wind in 
Fig. 2, except the 7-year quasi-periodicity. As stated earlier, these features presumably are part of the 
regional climatology associated with the North Atlantic Oscillation system. Note that the 6.5 month 
is close to half the 14-month periodicity and the 7-month quasi-periodicity reported by Naito & 
Kikuchi (1992) in global atmospheric angular momentum. 

The 7-year quasi-periodicity is relatively weak in our computed output, and certainly much 
weaker than the observed (see Fig. 9). The corresponding band, in fact, is the only frequency band 
where exists a major discrepancy in the spectral comparison of Fig. 9. The interesting feature, 
however, is that their relative behavior around this period seems highly correlated as revealed by the 
high coherence (Fig. 9) and similar wavelet spectra (Fig. 10). Thus, the causes for the discrepancy 
may be the excessive low-frequency smoothing of the NCEP stress field, or the existence of other 
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unmodeled forcing mechanism which is highly correlated with the wind-forcing such as air-sea 
temperature differences. In any event, this is outside of our present interest. A relatively smaller 
discrepancy in Fig. 9 exists around the 9-month peak. The above causes may play a role, but note that 
the relative size of the difference in amplitude is only on the order of 30%. At all other frequencies 
the computed sea level matches the observed quite well. This is to be further discuss below. 

5. CONCLUSIONS 

We have used the 40-year (1958-1997), monthly surface stress field in the North Sea region derived 
by NCEP as input to a barotropic North Sea model with realistic geography, bathymetry, and 
boundary conditions. The output wind-driven sea level variations are compared with 10 tide gauge 
records along the continental coast. Besides the strong seasonal variations, several prominent quasi- 
periodicities exist at around 7 years, 3 years, 14 months, 9 months, and 6.5 months. 

Correlation and spectral analyses show remarkable agreement between the model output and the 
observations, across almost the entire frequency band (out to the Nyquist frequency), except at the 
seasonal periods and periods longer than about 7 years. Our study only pertains to wind-driven sea 
level variations. The seasonal discrepancy presumably originates from other meteorological causes 
than wind-forcing that exist in the observation but are not modeled in our study. The long-period 
discrepancy may also have meteorological causes. 

Our focus here is the 14-month tide, which has been observed to be as large as 3 cm in amplitude 
along the Denmark-Netherlands coast. In the past both wind-forcing and certain North-Sea 
resonance-enhancement of the Chandler pole tide have been examined to explain this tide. Our 
present findings strongly supports the wind-forcing mechanism as the major cause of the 14-month 
tide. To be sure, a 14-month pole tide (and an annual pole tide for that matter) should in general be 
present; it’s equilibrium amplitude is no more than 0.5 cm. It was stated earlier that the computed sea 
level amplitude depends on the model parameterization and the wind stress forcing. The fact that in 
our simulation the amplitude matches across the entire frequency band indicates that the 
parameterization adopted here is realistic, and having no preference or “favor” towards the matching 
at the 14-month periodicity of interest here. In other words, we find no evidence to support an 
anomalous behavior of the North Sea in this particular band, and hence no need to invoke a 
geophysical explanation involving resonance-enhancement of pole tide to explain the observations. 
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Table 1. The ten North Sea tide gauge stations used in this study. The last two columns respectively 
list the correlation coefficients between the computed and the observed non-seasonal sea level signals, 
and the ratio [standard deviation of computed amplitude]/[standard deviation of observed amplitude]. 


Station 

Latitude 

Longitude 

Corr. Coef. 

STD 

ratio 

Esbjerg, Denmark 

55°N 28’ 

8°E 26' 

0.880 

0.90 

West-Terschelling, Netherlands 

53°N 22' 

5°E 13' 

0.784 

1.07 

Delfzijl, Netherlands 

53 °N 20' 

6°E 56' 

0.885 

0.95 

Harlingen, Netherlands 

53°N 10' 

5°E 25' 

0.833 

0.93 

Den Helder, Netherlands 

52°N 58' 

4°E 45' 

0.759 

1.11 

Ijmuiden, Netherlands 

52°N 28’ 

4°E 35' 

0.735 

1.26 

Hoek Van Holland, Netherlands 

51°N 59' 

4°E 7' 

0.730 

1.35 

Maassluis, Netherlands 

51°N 55' 

4°E 15' 

0.712 

1.24 

Vlissingen, Netherlands 

51 °N 27' 

3°E 36' 

0.674 

1.70 

Oostende, Belgium 

51°N 14' 

2°E 55' 

0.462 

1.76 
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Figure Captions: 

Figure 1. (a) Map of Europe showing the North Sea region modeled in this study, (b) The modeled 
North Sea region with depth contours in meters. The crosses indicate locations of the ten Permanent 
Service for Mean Sea Level tide gauge stations on the continental coast examined in this study (see 
Table 1), including the northern five Denmark-Netherlands coastal stations (marked by a), where the 
14-month tide is the strongest (and used to form the composite record later). Overlain as thick arrows 
are the 40-year (1958-1997) averaged wind stress field at NCEP grid points. 

Figure 2. (a) The 40- year time series of the horizontal wind velocity in the North Sea region, spatially 
averaged over the NCEP grid points; u: eastward component; v: northward component (shifted by 
-20 m/s in the plot), (b) Their maximum-entropy power spectra. 

Fig. 3. Sea level contours in centimeter, and currents driven by the 40-year average NCEP wind stress 
(Fig. 1), predicted by our North Sea ocean model. 

Figure 4. (a) Comparison of the annual signal of computed (*) wind-driven sea level with that 
observed (o), in terms of amplitude and phase as a function of the latitude of the ten coastal stations, 
(b) Same as (a) but for the semi-annual signals. 

Figure 5. Comparison of computed (top curves, vertically shifted) wind-driven sea level with that 
observed (bottom curves) after the removal of seasonal (annual and semi-annual) signals and a linear 
trend, shown for four randomly selected stations. 

Figure 6. Comparison of the computed (top curve, vertically shifted) and observed (bottom curve) 
“composite” average sea level of the five Denmark-Netherlands coastal stations, after the removal 
of seasonal (annual and semi-annual) signals and a linear trend. 

Figure 7. The (time-domain) cross-correlation function between the two series in Fig. 6. 

Figure 8. The (frequency-domain) coherence spectrum between the two series in Fig. 6. Upper panel: 
coherence amplitude squared compared with the 99% confidence level. Lower panel, coherence 
phase. 

Figure 9. The (multi- tapered) power spectra for the two series in Fig. 6. 

Fig. 10. The (time-frequency) wavelet spectra for the two series in Fig. 6 (upper panel: computed; 
lower panel observed). The contours indicate the strength of undulations at any particular frequency 
(or period) as a function of time. The brighter contour closures indicate positive amplitude and the 
darker ones the negative amplitude. 
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